body .main-container {
  max-width: 100% !important;
  width: 100% !important;
}
body {
  max-width: 100% !important;
}

body, td {
   font-size: 16px;
}
code.r{
  font-size: 14px;
}
pre {
  font-size: 14px
}
knitr::opts_chunk$set(message = FALSE, 
                      warning = FALSE, 
                      echo = TRUE, 
                      include = TRUE, 
                      fig.align = "center",
                      out.width = "100%"
                      )
set.seed(1982)
library(INLA)
library(here)

1 Model 1

data(Seeds)
df <- data.frame(y = Seeds$r, Ntrials = Seeds$n, Seeds[, 3:5])
df |> head(2) |> knitr::kable(caption = "Data")
Data
y Ntrials x1 x2 plate
10 39 0 0 1
23 62 0 0 2

1.1 Model specification

The model specification is

\[ y_i\sim\text{Binomial}\left(N_i, p_i\right)\label{eq1}\tag{A} \]

\[ p_i = \text{invlogit}(\eta_i) = \frac{\exp(\eta_i)}{1+\exp(\eta_i)} = \frac{1}{1+\exp(-\eta_i)}\qquad \left(\text{i.e., }\eta_i = \text{logit}(p_i) =\log\left(\frac{p_i}{1-p_i}\right) \right) \]

\[ \eta_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} +v_i \]

\[ \beta_i \sim \text{N}(0, 1000) \]

\[ v_i\sim\text{N}(0,\sigma^2) \] \[ \sigma\sim\text{Exp}(\lambda) \qquad \left(\text{i.e., }\pi(\sigma) = \lambda e^{-\lambda\sigma}\right) \]

\[ \texttt{inla()}\text{ deals with }\sigma^{-2} = \tau \sim\pi(\tau) =\pi(\sigma) \left|\dfrac{d\sigma}{d\tau}\right|= \dfrac{\lambda}{2}\tau^{-3/2}e^{-\lambda \tau^{-1/2}} \]

\[ \text{internally }\texttt{inla()}\text{ deals with }\log(\tau) = \theta \sim\pi(\theta) = \pi(\tau)\left|\dfrac{d\tau}{d\theta}\right| = \dfrac{\lambda}{2}\exp\left(-\dfrac{\theta}{2} - \lambda \exp\left(-\dfrac{\theta}{2}\right)\right) \]

\[ \lambda \text{ is such that } \pi(\sigma>1) = 0.01\qquad\left(\text{i.e., }0.01 = \pi(\sigma>1) = 1- \pi(\sigma\leq1) = e^{-\lambda}\implies \lambda = -\dfrac{\log(0.01)}{1}\approx4.6\right)\label{eq2}\tag{B} \]

1.2 INLA call

# call to inla
hyper <- list(theta = list(prior = "pc.prec", param = c(1,0.01)))
formula1 = y ~ x1 + x2 + f(plate, 
                           model = "iid",  
                           hyper = hyper)
model1 = inla(formula = formula1, 
            data = df, 
            family = "binomial",
            Ntrials = Ntrials, 
            control.family = list(control.link = list(model = "logit")),
            control.predictor = list(compute = TRUE, link = 1)) 
summary(model1)
## Time used:
##     Pre = 0.373, Running = 0.159, Post = 0.0129, Total = 0.545 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.387 0.181     -0.739   -0.390     -0.017 -0.390   0
## x1          -0.360 0.231     -0.838   -0.353      0.076 -0.353   0
## x2           1.033 0.221      0.589    1.034      1.469  1.034   0
## 
## Random effects:
##   Name     Model
##     plate IID model
## 
## Model hyperparameters:
##                      mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for plate 19.78 32.61       2.96    10.55      84.94 6.06
## 
## Marginal log-Likelihood:  -68.44 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

2 Model 2

# data generation
N = 3600

sig.epsilon = 0.5
sig.u = 1
sig.seasonal = 2

rho = 0.90

u = arima.sim(list(order = c(1, 0 ,0), ar = rho), n = N,sd = 1)
u = u/sd(u)*sig.u
seas.coeff = (0:11)*(1:12-12)
seas = rep(seas.coeff, N/12)
seas = drop(scale(seas))*sig.seasonal
gaussian.error = rnorm(N, 0, sd = sig.epsilon)
(prec.sig.epsilon = sig.epsilon^(-2))
## [1] 4
(prec.sig.u = sig.u^(-2))
## [1] 1
(prec.sig.seasonal = sig.seasonal^(-2))
## [1] 0.25
# the model
y = u + seas + gaussian.error
df = data.frame(y = y, t = 1:N, year = rep(1:12, N/12))
df |> head(2) |> knitr::kable(caption = "Data")
Data
y t year
1.9358080 1 1
0.8643198 2 2

2.1 Model specification

\[ y_i = \beta_0 + u(t_i) + s_i + \epsilon_i,\; i =1,\dots, n \]

\[ u(t_1) \sim \text{N}(0, (\tau_u(1-\rho^2))^{-1}), \qquad u(t_i) = \rho u(t_{i-1}) + e_i, \qquad e_i \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \qquad i =2,\dots, n, \qquad |\rho|<1\qquad \text{or}\qquad u(t_i) - \rho u(t_{i-1}) \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \]

\[ \beta_0 \sim \text{N}(0, 1000) \]

\[ \epsilon_i \sim \text{N}(0, \tau_\epsilon^{-1} = \sigma_\epsilon^2) \]

\[ s_i-2s_{i+1}+s_{i+2}\sim\text{N}(0, \tau_s^{-1} = \sigma_s^2) \]

\[ \sigma_u^{-2} = \tau_u \sim\pi(\tau_u) = \dfrac{\lambda}{2}\tau_u^{-3/2}e^{-\lambda \tau_u^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.1}\text{ is such that }\pi(\sigma_u>0.1) = 0.5 \]

\[ \rho\sim\pi(\rho) = \dfrac{\lambda\exp(-\lambda\sqrt{1-\rho})}{1-\exp(-\sqrt2\lambda)}\dfrac{1}{2\sqrt{1-\rho}},\qquad\text{where }\lambda,\text{ the solution of }\dfrac{\exp(-\lambda\sqrt{1-0.9})}{1-\exp(-\sqrt2\lambda)} = 0.5,\text{ is such that }\pi(\rho>0.9) = 0.5 \]

\[ \sigma_s^{-2} = \tau_s \sim\pi(\tau_s) = \dfrac{\lambda}{2}\tau_s^{-3/2}e^{-\lambda \tau_s^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.1}\text{ is such that }\pi(\sigma_s>0.1) = 0.5 \]

\[ \sigma_\epsilon^{-2} = \tau_\epsilon \sim\pi(\tau_\epsilon) = \dfrac{\lambda}{2}\tau_\epsilon^{-3/2}e^{-\lambda \tau_\epsilon^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{3}\text{ is such that }\pi(\sigma_\epsilon>3) = 0.5 \]

2.2 INLA call

hyper.ar1 <- list(theta1 = list(prior="pc.prec", param = c(0.1, 0.5)), #sigma_u
                 theta2 = list(prior="pc.cor1", param = c(0.9, 0.5))) #rho
hyper.rw2 <- list(theta1 = list(prior="pc.prec", param = c(1.99, 0.99))) #sigma_s
hyper.family <- list(theta = list(prior="pc.prec", param = c(3, 0.5))) #sigma_epsilon

formula <- y ~ f(t, model = 'ar1', hyper = hyper.ar1) + 
  f(year, model = "rw2", hyper = hyper.rw2, cyclic = T, constr = T)

model2 <- inla(formula = formula,
            data = df, 
            family = "gaussian",
            control.predictor = list(compute = TRUE),
            control.family = list(hyper = hyper.family))
summary(model2)
## Time used:
##     Pre = 0.192, Running = 1.1, Post = 0.164, Total = 1.45 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.043 0.076     -0.106    0.043      0.192 0.043   0
## 
## Random effects:
##   Name     Model
##     t AR1 model
##    year RW2 model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.990 0.184      3.640    3.986
## Precision for t                         1.018 0.078      0.870    1.016
## Rho for t                               0.907 0.009      0.889    0.908
## Precision for year                      1.210 0.546      0.439    1.111
##                                         0.975quant  mode
## Precision for the Gaussian observations      4.365 3.977
## Precision for t                              1.178 1.013
## Rho for t                                    0.924 0.908
## Precision for year                           2.542 0.929
## 
## Marginal log-Likelihood:  -4048.59 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
formula <- y ~ f(t, model = 'ar1') + 
  f(year, model = "rw2", cyclic = T, constr = T)

model2 <- inla(formula = formula,
            data = df, 
            family = "gaussian",
            control.predictor = list(compute = TRUE))
summary(model2)
## Time used:
##     Pre = 0.215, Running = 1.49, Post = 0.274, Total = 1.98 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.043 0.076     -0.107    0.043      0.194 0.043   0
## 
## Random effects:
##   Name     Model
##     t AR1 model
##    year RW2 model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.997 0.185      3.642    3.994
## Precision for t                         1.007 0.078      0.859    1.005
## Rho for t                               0.908 0.009      0.890    0.908
## Precision for year                      1.570 0.633      0.651    1.464
##                                         0.975quant  mode
## Precision for the Gaussian observations      4.371 3.990
## Precision for t                              1.167 1.003
## Rho for t                                    0.925 0.908
## Precision for year                           3.099 1.270
## 
## Marginal log-Likelihood:  -4057.16 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

3 Model 3

# temp = read.csv("data/harmonised-unemployment-rates-mo.csv")
n = nrow(temp) - 1
df = data.frame(y = temp[1:n, 2], t = 1:n, dates = temp[1:n, 1])
df |> head(2) |> knitr::kable(caption = "Data")
Data
y t dates
9.1 1 2000-03
8.8 2 2000-04

3.1 Model specification

\[ y_i = \beta_0+u(t_i)+\epsilon_i \]

\[ u(t_1) \sim \text{N}(0, (\tau_u(1-\rho^2))^{-1}), \qquad u(t_i) = \rho u(t_{i-1}) + e_i, \qquad e_i \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \qquad i =2,\dots, n, \qquad |\rho|<1\qquad\text{or}\qquad u(t_i) - \rho u(t_{i-1}) \sim \text{N}(0, \tau_u^{-1} = \sigma_u^2) \]

\[ \beta_0 \sim \text{N}(0, 1000) \]

\[ \epsilon_i \sim \text{N}(0, \tau_\epsilon^{-1} = \sigma_\epsilon^2) \]

\[ \sigma_u^{-2} = \tau_u \sim\pi(\tau_u) = \dfrac{\lambda}{2}\tau_u^{-3/2}e^{-\lambda \tau_u^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{0.02}\text{ is such that }\pi(\sigma_u>0.02) = 0.5 \]

\[ \rho\sim\pi(\rho) = \dfrac{\lambda\exp(-\lambda\sqrt{1-\rho})}{1-\exp(-\sqrt2\lambda)}\dfrac{1}{2\sqrt{1-\rho}},\qquad\text{where }\lambda,\text{ the solution of }\dfrac{\exp(-\lambda\sqrt{1-0.9})}{1-\exp(-\sqrt2\lambda)} = 0.5,\text{ is such that }\pi(\rho>0.9) = 0.5 \]

\[ \sigma_\epsilon^{-2} = \tau_\epsilon \sim\pi(\tau_\epsilon) = \dfrac{\lambda}{2}\tau_\epsilon^{-3/2}e^{-\lambda \tau_\epsilon^{-1/2}},\qquad\text{where }\lambda = -\dfrac{\log(0.5)}{3}\text{ is such that }\pi(\sigma_\epsilon>3) = 0.5 \]

Notice that \(\sigma_u\) and \(\sigma_\epsilon\) are \(\text{Exp}(\lambda)\).

3.2 INLA call

hyper.ar1 = list(theta1 = list(prior = "pc.prec", param = c(0.02, 0.5)), #sigma_u
                 theta2 = list(prior = "pc.cor1", param = c(0.9, 0.5))) #rho
hyper.family = list(theta = list(prior = "pc.prec", param = c(3, 0.5))) #sigma_epsilon
formula2 <- y ~ f(t,model = 'ar1', hyper = hyper.ar1)

model3 <- inla(formula = formula2,
             data = df,
             family="gaussian",
             control.predictor = list(compute=TRUE),
             control.family = list(hyper = hyper.family))
summary(model3)
## Time used:
##     Pre = 0.149, Running = 0.194, Post = 0.0109, Total = 0.354 
## Fixed effects:
##             mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 8.64 0.305      8.034    8.642      9.239 8.642   0
## 
## Random effects:
##   Name     Model
##     t AR1 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.20e+04 2.42e+04   1452.689 1.44e+04
## Precision for t                         6.12e-01 7.70e-02      0.472 6.09e-01
## Rho for t                               7.84e-01 2.80e-02      0.726 7.86e-01
##                                         0.975quant     mode
## Precision for the Gaussian observations   8.64e+04 3952.496
## Precision for t                           7.76e-01    0.603
## Rho for t                                 8.36e-01    0.788
## 
## Marginal log-Likelihood:  -236.96 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

4 Model 4